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Abstract Single- and multi-band Hubbard models have 
been found to describe many of the complex phenomena 
that are observed in the cuprate and iron-based high- 
temperature superconductors. Simulations of these mod- 
els therefore provide an ideal framework to study and 
understand the superconducting properties of these sys- 
tems and the mechanisms responsible for them. Here we 
review recent dynamic cluster quantum Monte Carlo 
simulations of these models, which provide an unbiased 
view of the leading correlations in the system. In par- 
ticular, we discuss what these simulations tell us about 
superconductivity in the homogeneous 2D single-orbital 
Hubbard model, and how charge stripes affect this be- 
havior. We then describe recent simulations of a bi- 
layer Hubbard model, which provides a simple model 
to study the type and nature of pairing in systems with 
multiple Fermi surfaces such as the iron-based super- 
conductors. 
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Phenomenologically, the cuprate and iron-based high- 
temperature superconductors share in common that su- 
perconductivity occurs by doping into an antiferromag- 
netic Mott or spin density wave state. The magnetism 
in these compounds originates from partially filled Cu 
or Fe d-orbitals which energetically lay at the Fermi 
energy. A reasonable starting point for a theoretical 
description of these compounds is therefore given by 
a Hubbard model which describes moment formation 
due to a strong local Coulomb repulsion between the 
electrons on the d-orbitals. In the case of the cuprates, 
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a single electronic band mainly of Cu d^^^-character 
crosses the Fermi surface and a single-band Hubbard 
model has therefore been argued to provide a simple 
framework to describe the low-energy physics of this 
class of materials [Til 1 8 J . In the iron-based materials, 
one has multiple Fermi surfaces formed by Bloch states 
originating from several of the iron d-orbitals. There- 
fore, one needs to start with a multi-orbital Hubbard 
model to describe the more complex electronic struc- 
ture of these systems [51ITU] . A two-band model has also 
been used to study phase separation in cuprate super- 
conductors [5]. 

Here, we review recent dynamic cluster quantum 
Monte Carlo calculations of these models [14l[T5lfr2l[T6] . 
In particular, we will focus on what these simulations 
tell us about the nature of pairing and the effect of 
charge stripes in the single-band 2D Hubbard model, 
as well as the leading pairing correlations in a bilayer 
Hubbard model with multiple Fermi surfaces. 

The Hamiltonian for the 2D Hubbard model we will 
study can be written as 



n = -tJ2 (4 



C 

la ja 



h.c.) + Uj2 n 



(1) 



Here t is a nearest neighbor hopping parameter, U is 
an on-site Coulomb repulsion and (ij) implies summa- 
tion over nearest neighbor pairs only. In the following 
we will measure energies in units of t. The single band 
with dispersion Ek = —2t(cosk x +cosk v ) describes the 
electronic states near the Fermi energy in the cuprate 
materials. 

As discussed in the introduction, a realistic descrip- 
tion of the pnictides requires a multi-orbital Hubbard 
model describing the 5 iron d-orbitals near the Fermi 
energy and their intra- and inter-orbital Coulomb, Hund's 
rule and pair hopping interactions. Such a model is too 
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complex to simulate with quantum Monte Carlo ap- 
proaches, and so we will instead study a simpler two- 
orbital model with only intra-orbital Coulomb interac- 
tions. This model will be realized by a bilayer Hubbard 
model. Its Hamiltonian is given by 
\ " / t „ 



H = -t 



(ij)ma 



^"im'j'^'imj, ■ 



+ hxO-t±y>Lc i2(r + h.c.) 



(2) 



Here the layers are indexed by m, each layer is de- 
scribed by the Hamiltonian in Eq. [T] and t± is an ad- 
ditional hopping parameter between neighboring sites 
in the bi-layer model. This model provides a simplified 
two-orbital system in which one can study the type of 
pairing that can occur in systems with multiple Fermi 
surfaces such as the iron-pnictides. 

In order to analyze these models, we will use a dy- 
namic cluster quantum Monte Carlo approximation 
(DCA/QMC) 0SCE3]. The DCA maps the bulk lattice 
problem onto an effective periodic cluster cluster em- 
bedded in a dynamic mean-field that is designed to rep- 
resent the rest of the system. The effective cluster prob- 
lem is then solved using a quantum Monte Carlo algo- 
rithm. The results discussed in this paper were obtained 
with a Hirsch-Fye method [5]. DCA/QMC calculations 
of the 2D Hubbard model have found many phenom- 
ena that are also observed in the cuprates, including 
an antiferromagnetic Mott state, d-wave superconduc- 
tivity as well as pseudogap behavior [T3]. It therefore 
provides an interesting framework to study many of the 
open questions in the field. 

Formally, the quantity of interest to study the na- 
ture of pairing in these models is given by the two- 
particle irreducible vertex in the particle-particle chan- 
nel, r?£(k, k') |14) . Here k = (k, ioj n ) with w„ a fermion 
Matsubara frequency and we are interested in the sin- 
glet pairing channel. This quantity describes the scat- 
tering of two electrons with momenta k and — k and 
anti-parallel spins to a state with momenta k' and — k' 
and therefore describes the pairing interaction. Together 
with the single-particle Green's function G(k), it enters 
the Bethe-Salpeter equation 



^E^(M')G(fc')G(-fc')<Mfc') 

k 



XMk) (3) 



which provides information on the strength (A Q ) and 
momentum and frequency structure (^(fc)) of the lead- 
ing pairing correlations in the system [H]. At T c , X a = 1 
and ^q,(A;) becomes identical to the superconducting 
gap. In the 2D Hubbard model, at low temperature, one 
finds that the eigenvector corresponding to the leading 
eigenvalue has a d-wave cos k x 
pendence. 



cos k y momentum de- 



Prcvious DCA/QMC simulations of the 2D Hub- 
bard model [H¥T5] have found that the momentum 
and frequency dependence of the pairing interaction 
r?P(k,k') is similar to that of the spin susceptibility 
x{k — k'), providing evidence that that pairing interac- 
tion in this model is carried by spin fluctuations. 

In a spin fluctuation picture, one can naturally un- 
derstand the drop of T c with doping on the overdoped 
side of the cuprate phase diagram, since the spin-fluctuations 
are weakened by doping away from the antiferromag- 
netic parent state. On the other hand, the drop of T c 
with underdoping is difficult to understand in a picture 
where pairing is mediated by spin fluctuations, since 
one would expect them to get stronger when the sys- 
tem is doped towards the Mott state. To investigate 
this issue we show in Fig. la the temperature versus 
doping superconducting phase diagram of the 2D Hub- 
bard model, calculated with DCA/QMC on an 8-site 
cluster with U=8. One sees that these calculations cor- 
rectly predict the experimentally observed dome-like 
structure of the superconducting phase diagram, with 
T c dropping with both over- and underdoping. 
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Fig. 1 (a) Superconducting transition temperature T c ver- 
sus doping x for the 2D Hubbard model with (7 = 8 calculated 
with DCA/QMC on an 8-site cluster, (b) Normalized interac- 
tion strength Vd and "intrinsic" pair- field susceptibility Pd,o 
versus doping x calculated at a temperature T = 0.125. 



In order to analyze how this behavior arises from 
the Bethe-Salpeter equation ©, we have calculated the 
"intrinsic" pair- field susceptibility projected onto the 
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leading eigenvector, P dfi = T/N J2 k $d{k) 2 G(k)G(-k) 
and the strength of the pairing interaction Vd from 
VdPd,o = ^d- The doping dependence of these quanti- 
ties calculated at a low temperature above T c are shown 
in Fig. lb. As one would expect from a spin- fluctuation 
based pairing interaction, Vd rises monotonically with 
decreasing doping towards the Mott insulator. In con- 
trast, Pd : o decreases with decreasing doping and goes to 
zero as one approaches the Mott state. The drop of T c 
on the underdoped side therefore is caused by the strong 
Mott quasiparticle renormalization when the doping is 
close to zero. 

Despite the strong pairing interaction Vd in the un- 
derdoped system, T c is small because of the Mott quasi- 
particle degradation. It has been suggested that in this 
case, T c may be enhanced in a striped state, in which 
the system is separated into hole-rich regions with good 
hole mobility, and hole-poor regions with strong anti- 
ferromagnetic correlations. To explore this prediction, 
we have performed DCA/QMC calculations of a striped 
8x4-site cluster [TJ]. The charge stripes were imposed 
by hand, by adding a site-dependent external potential 
Vi that couples to the charge density on site i. Vi is 
chosen to have a maximum amplitude Vq and to vary 
smoothly along the long 8-site x-direction of the cluster 
while being constant along the 4-sitc y-dircction. The 
variation of Vi along the x-direction is shown at the 
bottom of Fig. 2a. Here we are interested in a striped 
inhomogeneity with period 8. The calculated site- filling 
is displayed at the top of Fig. 2a and found to follow 
closely the variation of the external potential. 

In order to keep the problem computationally tractable, 
we have averaged over different stripe origins along the 
x-direction, before computing the mean-field medium 
[TJ] . This corresponds to a situation where the stripe or- 
der is short-ranged, over the length-scale of the cluster, 
but is translationally invariant on longer macroscopic 
length scales. 

The temperature dependence of the leading (d-wave) 
eigenvalue of the Bethe-Salpeter equation (J2|) computed 
for different amplitudes Vq of the potential is shown in 
Fig. 2b. As one can see, the pairing correlations are 
indeed enhanced by the charge stripe and T c , i.e. the 
temperature where Xd(T) crosses 1, is increased. One 
also sees that there is an optimum inhomogeneity at 
around Vq — 0.4 for which T c is maximized. For larger 
Vq, T c is found to drop again. 

Finally, we discuss recent DCA/QMC results for the 
bilayer Hubbard model [16] . This model provides a sim- 
ple system with multiple Fermi surfaces, in which one 
can study the type of pairing that can occur in sys- 
tems such as the iron-pnictides. Due to the inter-layer 
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Fig. 2 DCA/QMC simulation of an inhomogeneous 8x4 
cluster with a charge stripe, (a) External potential V(l x ) and 
resulting charge density (n)i x along the long (8-site) direc- 
tion of the cluster for different magnitudes Vq. (b) Leading 
eigenvalue of the particle-particle Bethe-Salpeter equation ([3J 
versus temperature for different Vq. 



hopping t±, one has a bonding and an anti-bonding 
band split by 2t± for U = 0. The bonding and anti- 
bonding Fermi surfaces of the non-interacting system 
are shown in Fig. 3 (a) and (b) for a filling (n) = 0.95 
and two values of t±/t. As t±/t increases the Fermi sur- 
faces shrink and for (n) = 1 the non-interacting system 
becomes a band-insulator for t± > 4. Weak coupling 
calculations for the doped system [3lHT] have found a 
(4,2 _ y 2 -like gap for t±/t = 0.5 and a fairly isotropic 
"s gap that changes sign between the bonding and 
anti-bonding Fermi surfaces for t±/t = 2.0, as schemat- 
ically illustrated in Fig. 3 (a) and (b). 

Here we describe the results of DCA/QMC calcula- 
tions which provide an unbiased treatment of this model 
and allow us to reach lower temperatures than previous 
finite size quantum Monte Carlo calculations [31H7U71H1 
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of the pairfield susceptibilities can be understood in 
terms of the £j_-dependence of the spin-fluctuation spec- 
tral weight. For small t±, the intra-layer spin fluctua- 
tions are dominant and give rise to the d x 2_ y 2 response. 
As t± increases the dominant spin-fluctuations change 
from intra- to inter-layer, which give rise to the s ± pair- 
ing. For large ij_ > 3, the inter-layer spin fluctuations 
become gapped and contribute less to the pairing. 

To conclude, we have reviewed dynamic cluster quan- 
tum Monte Carlo simulations of models of unconven- 
tional superconductors, including a 2D Hubbard model 
with charge stripes and a bilayer Hubbard model with 
multiple Fermi surfaces. We have shown that charge 
stripes in the 2D Hubbard model can lead to a sig- 
nificant enhancement of superconductivity. For the bi- 
layer model, we have found a transition of the leading 
pairing instability from a d x 2_ y 2- wave to an s ± state 
with increasing inter-layer hopping t±/t. We have also 
discussed how the superconducting behavior of these 
models can be understood in terms of a spin-fluctuation 
picture. 



Fig. 3 DCA/QMC simulation of a bilayer Hubbard model 
with inter-layer hopping t±. The bonding (solid lines) and 
anti-bonding (dashed lines) Fermi surfaces of the non- 
interacting system for (a) t± = 0.5 and (b) t± = 2 for a 
filling (n) = 0.95. A d x 2_ y 2 gap structure is illustrated in 
for the t±/t — 0.5 Fermi surface and an gap structure is 
shown for t±/t = 2.0. (c) The d-wave and pair-field sus- 
ceptibilities P a versus temperature for different values of the 
inter-layer hopping t±. As t± increases, the leading instabil- 
ity changes from d-wave to and at larger values of t± the 
superconducting pair-field susceptibility is suppressed. 
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12]. The superconducting response is studied by calcu- 
lating the pairfield susceptibility 

/>oo 

P a (T)= / dr(A a (r)AU0)), (4) 
Jo 

Here A^ = 3(^) c lt c -fe4. an< ^ = cos ^x — 

cos k y for the d x 2_ y 2-w&ve and cosfc z for the s ± case. 
In Fig. 3 (c)-(f), results for these pair- field suscepti- 
bilities versus temperature are shown for different val- 
ues of t±. For t± = 0.5, the d-wave response is larger 
than the s ± susceptibility. For t± = 1, both channels 
are almost degenerate, while for t± = 2, the s is the 
leading response and diverges at a relatively high tem- 
perature around 0.05t. For tj_ = 3, the pair-field re- 
sponse is significantly weakened. This near-degeneracy 
of the d x 2_ y 2 and pair-field correlations for inter- 
mediate values of t±/t has been similarly observed in 
fluctuation-exchange calculations of realistic 5-orbital 
model calculations of the iron-pnictides [SlfTU] . 

A further analysis of the pairing interaction in the 
bilayer model [T^ shows that the observed ^-dependence 
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